## Visualize experiment from hives 7-8 ## Bees were given two treatments -- reward first or second
#install packages
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp",
"plyr", "effects", "viridis", "GGally", "survminer",
"tidyr","gamm4")
ipak(packages)## Loading required package: ggplot2## Loading required package: reshape2## Loading required package: lme4## Loading required package: Matrix## Loading required package: sjPlot## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!## Loading required package: multcomp## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: plyr## Loading required package: effects## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Loading required package: viridis## Loading required package: viridisLite## Loading required package: GGally## Loading required package: survminer## Loading required package: ggpubr## Loading required package: magrittr## Loading required package: tidyr##
## Attaching package: 'tidyr'## The following object is masked from 'package:magrittr':
##
## extract## The following object is masked from 'package:Matrix':
##
## expand## The following object is masked from 'package:reshape2':
##
## smiths## Loading required package: gamm4## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.## This is gamm4 0.2-5## ggplot2 reshape2 lme4 sjPlot multcomp plyr effects
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## viridis GGally survminer tidyr gamm4
## TRUE TRUE TRUE TRUE TRUE# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))
# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"
# print system info
print(paste("last run ", Sys.time()))## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'## [1] "last run 2017-12-09 06:53:18"print(R.version)## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 4.2
## year 2017
## month 09
## day 28
## svn rev 73368
## language R
## version.string R version 3.4.2 (2017-09-28)
## nickname Short SummerRead in data
bees <- read.csv(file.path(dataDir, "03_StartStop_cleaned.csv"))
bees$hive = as.factor(bees$hive)
head(bees)## freq amp datetime visitNum treatment_rewarded
## 1 370 0.79858 2017_05_30__09_52_48_207 1 F
## 2 340 0.68663 2017_05_30__09_52_50_473 2 F
## 3 330 0.34916 2017_05_30__09_52_52_817 3 F
## 4 320 0.66814 2017_05_30__09_52_56_614 4 F
## 5 310 0.90015 2017_05_30__09_52_58_795 5 F
## 6 320 1.10685 2017_05_30__09_52_59_487 6 F
## BeeNumCol accFile
## 1 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_48_207_220_450_0_5.txt
## 2 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_50_473_220_450_0_5.txt
## 3 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_52_817_220_450_0_5.txt
## 4 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_56_614_220_450_0_5.txt
## 5 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_58_795_220_450_0_5.txt
## 6 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_59_487_220_450_0_5.txt
## beeID hive treatment totalRewards meanFreqRewarded
## 1 207 7 RewSec 50 303.8
## 2 207 7 RewSec 50 303.8
## 3 207 7 RewSec 50 303.8
## 4 207 7 RewSec 50 303.8
## 5 207 7 RewSec 50 303.8
## 6 207 7 RewSec 50 303.8
## totalBuzzesNotRewarded meanFreqUnrewarded IT_mm amp_acc
## 1 50 307.2 4.11 78.52311
## 2 50 307.2 4.11 67.51524
## 3 50 307.2 4.11 34.33235
## 4 50 307.2 4.11 65.69715
## 5 50 307.2 4.11 88.51032
## 6 50 307.2 4.11 108.83481# calculate mean freq rewarded vs. unrewarded for each bee
# to make sure it agrees with what is already in the dataframe
tapply(bees$freq, INDEX = list(bees$beeID, bees$treatment_rewarded), mean)## F T
## 207 307.2000 303.8000
## 208 344.6000 338.0000
## 209 317.8000 334.6000
## 210 342.5000 343.6000
## 211 376.3265 NA
## 212 331.5152 NA
## 213 NA 337.2093
## 214 327.8000 349.8000
## 215 NA 356.6667
## 216 366.8000 342.6000
## 217 310.0000 NA
## 218 363.2000 327.0000
## 219 353.6000 312.6000
## 220 340.4000 336.4000
## 221 358.1818 NA
## 222 331.0000 304.4000
## 223 347.5000 NA
## 224 350.8000 348.8000
## 225 347.7419 354.4000
## 226 NA 361.4634
## 227 307.4000 309.3878
## 228 336.2000 361.4000
## 229 361.2000 351.0000
## 230 336.2000 298.2000
## 231 347.8000 367.4000
## 232 358.0000 343.8000
## 233 363.5714 NA
## 234 280.0000 313.0000
## 235 NA 352.0588
## 236 394.4828 NA
## 237 333.2000 313.2000
## 238 305.3333 NA
## 239 401.4000 395.0000
## 240 342.3810 NA
## 241 343.2000 319.8000
## 242 303.1579 NA
## 243 346.2000 348.0000
## 244 380.0000 352.2000
## 245 296.3636 NA
## 246 352.4000 329.1429
## 247 314.0000 304.6000
## 248 389.2000 349.8000
## 249 346.8000 336.4000
## 250 355.6000 347.4000
## 251 364.6341 NA
## 252 299.3333 NA
## 253 380.0000 NA
## 254 321.3636 314.0000
## 255 347.2000 NA
## 256 273.7500 NA
## 257 320.0000 347.8000
## 258 302.2222 NA
## 259 330.0000 NA
## 260 NA 328.2143
## 261 328.2353 NA
## 262 343.3333 NA
## 263 336.8000 315.2000
## 264 NA 294.0000
## 265 NA 311.3043
## 266 NA 356.6667
## 267 334.5714 330.4000
## 268 328.9362 NA
## 269 329.6000 363.6000
## 270 339.6000 344.2000
## 271 333.2609 NA
## 272 327.6471 NA
## 273 357.6471 357.4000
## 274 275.5556 NA
## 275 345.6000 336.6000
## 276 NA 336.4286
## 277 353.8000 333.2000
## 278 319.0000 305.2000
## 279 307.6000 287.8000
## 280 295.0000 NA
## 281 NA 344.2105
## 282 355.2000 360.0000
## 283 356.1905 NA
## 284 323.8462 NA
## 285 345.2000 340.8000
## 286 305.4000 297.2000
## 287 NA 386.2857
## 288 348.0000 345.1724
## 289 340.9091 NA
## 290 386.4000 359.0000
## 291 401.4706 369.8000
## 292 NA 347.7273
## 293 368.9474 NA
## 294 368.6000 365.8000
## 295 331.0526 NA
## 296 316.6000 312.4000
## 297 411.7778 NA
## 298 328.4000 314.2000
## 299 NA 318.8889
## 300 297.4000 289.8000
## 301 NA 385.8824
## 302 381.2903 NA# number of bees -- should be 96
length(unique(bees$BeeNumCol))## [1] 96# num observations -- should be 5851
nrow(bees)## [1] 5851Visualize data
ggplot(bees, aes(y = freq, x = treatment_rewarded)) +
geom_boxplot() +
stat_smooth(method = "loess", span = 0.75)
bees$trt2 <- mapvalues(bees$treatment, from = c("RewFir", "RewSec"),
to = c("Rewarded 1st 50 Buzzes", "Rewarded 2nd 50 Buzzes"))
ggplot(bees, aes(x = visitNum, y = freq, color = treatment_rewarded)) +
stat_smooth(method = "loess", span = 1) +
xlab("Buzz Number") +
ylab("Buzz Frequency(hz)") +
facet_wrap(~trt2) +
scale_color_viridis(name = "Rewarded Buzzes",
discrete = TRUE, begin = 0.2, end = 0.6)
#ggsave("beeRewards.png", width = 8, height = 5)
# make survival plot
survBees <- t(sapply(unique(bees$beeID), function(x) {
tmp <- bees[bees$beeID == x, ]
return(tmp[which.max(tmp$visitNum), ])
}))
lastObs <- as.data.frame(survBees)
lastObs$index <- as.numeric(lastObs$visitNum)
# 1 = censored, 2 = stopped buzzing for 5 min
lastObs$status <- mapvalues(lastObs$index == 100, from = c(TRUE, FALSE), to = c(1, 2))
# create survival object
lastObs$survObj <- with(lastObs, Surv(index, status == 2))
head(lastObs)## freq amp datetime visitNum treatment_rewarded BeeNumCol accFile
## 1 300 0.88263 100 100 2 1 100
## 2 320 0.90242 200 100 2 2 200
## 3 350 0.28458 300 100 1 3 300
## 4 310 0.84888 374 74 1 4 374
## 5 400 0.46176 423 49 1 5 423
## 6 330 0.57663 456 33 1 6 456
## beeID hive treatment totalRewards meanFreqRewarded
## 1 207 1 2 50 303.8
## 2 208 1 2 50 338
## 3 209 1 1 50 334.6
## 4 210 1 1 50 343.6
## 5 211 2 2 0 NA
## 6 212 2 2 0 NA
## totalBuzzesNotRewarded meanFreqUnrewarded IT_mm amp_acc trt2 index
## 1 50 307.2 4.11 86.78761 2 100
## 2 50 344.6 4.46 88.73353 2 100
## 3 50 317.8 4.4 27.9823 1 100
## 4 24 342.5 4.28 83.46903 1 74
## 5 49 376.3265 4.96 45.40413 2 49
## 6 33 331.5152 4.39 56.69912 2 33
## status survObj
## 1 1 100+
## 2 1 100+
## 3 1 100+
## 4 2 74
## 5 2 49
## 6 2 33lastObs$treatment <- as.factor(unlist(lastObs$treatment))
## Kaplan-Meier estimator.
km.as.one <- survfit(survObj ~ 1, data = lastObs, conf.type = "log-log")
km.by.trt <- survfit(survObj ~ treatment, data = lastObs, conf.type = "log-log")
par(mfrow = c(1,1))
plot(km.by.trt)
gsv <- ggsurvplot(km.by.trt, CI = FALSE)
gsv$plot + facet_wrap(~treatment) +
theme_bw() +
theme(legend.position = "none") +
scale_color_viridis(discrete = TRUE, begin = 0.3, end = 0.7, option = "magma") +
xlab("Buzz Number") +
ylab("Proportion Sonicating") 
#ggsave("beeSonicationPercentage.pdf", width = 8, height = 5)
ggsurv(km.by.trt, CI = FALSE, surv.col = c("#a6cee3", "#1f78b4")) +
ylim(c(0, 1)) +
xlab("Buzz Number") +
ylab("Proportion Sonicating") ## Loading required package: scales##
## Attaching package: 'scales'## The following object is masked from 'package:viridis':
##
## viridis_pal
xtabs(~lastObs$treatment + as.character(lastObs$hive))## as.character(lastObs$hive)
## lastObs$treatment 1 2
## RewFir 30 20
## RewSec 23 23Visualize data – averaged by bee
This is not super informative
beeAvgs = bees[bees$visitNum == 1, ]
nrow(beeAvgs) # should be 96## [1] 96head(beeAvgs)## freq amp datetime visitNum treatment_rewarded
## 1 370 0.79858 2017_05_30__09_52_48_207 1 F
## 101 350 1.25158 2017_05_30__10_15_55_623 1 F
## 201 390 1.19220 2017_05_30__10_28_20_624 1 T
## 301 350 0.79177 2017_05_30__10_38_37_411 1 T
## 375 390 0.78658 2017_05_30__11_13_02_616 1 F
## 424 270 0.04932 2017_05_30__11_27_40_579 1 F
## BeeNumCol accFile
## 1 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_48_207_220_450_0_5.txt
## 101 B_208_2017_05_30_H_7_RewSec 2017_05_30__10_15_55_623_220_450_0_5.txt
## 201 B_209_2017_05_30_H_7_RewFir 2017_05_30__10_28_20_624_220_450_0_5.txt
## 301 B_210_2017_05_30_H_7_RewFir 2017_05_30__10_38_37_411_220_450_0_5.txt
## 375 B_211_2017_05_30_H_8_RewSec 2017_05_30__11_13_02_616_220_450_0_5.txt
## 424 B_212_2017_05_30_H_8_RewSec 2017_05_30__11_27_40_579_220_450_0_5.txt
## beeID hive treatment totalRewards meanFreqRewarded
## 1 207 7 RewSec 50 303.8
## 101 208 7 RewSec 50 338.0
## 201 209 7 RewFir 50 334.6
## 301 210 7 RewFir 50 343.6
## 375 211 8 RewSec 0 NA
## 424 212 8 RewSec 0 NA
## totalBuzzesNotRewarded meanFreqUnrewarded IT_mm amp_acc
## 1 50 307.2000 4.11 78.523107
## 101 50 344.6000 4.46 123.065880
## 201 50 317.8000 4.40 117.227139
## 301 24 342.5000 4.28 77.853491
## 375 49 376.3265 4.96 77.343166
## 424 33 331.5152 4.39 4.849558
## trt2
## 1 Rewarded 2nd 50 Buzzes
## 101 Rewarded 2nd 50 Buzzes
## 201 Rewarded 1st 50 Buzzes
## 301 Rewarded 1st 50 Buzzes
## 375 Rewarded 2nd 50 Buzzes
## 424 Rewarded 2nd 50 Buzzesdata_long <- gather(beeAvgs[, c("BeeNumCol", "meanFreqRewarded", "meanFreqUnrewarded")], condition, meanFreq, meanFreqRewarded, meanFreqUnrewarded)
head(data_long)## BeeNumCol condition meanFreq
## 1 B_207_2017_05_30_H_7_RewSec meanFreqRewarded 303.8
## 2 B_208_2017_05_30_H_7_RewSec meanFreqRewarded 338.0
## 3 B_209_2017_05_30_H_7_RewFir meanFreqRewarded 334.6
## 4 B_210_2017_05_30_H_7_RewFir meanFreqRewarded 343.6
## 5 B_211_2017_05_30_H_8_RewSec meanFreqRewarded NA
## 6 B_212_2017_05_30_H_8_RewSec meanFreqRewarded NAggplot(data_long, aes(x = condition, y = meanFreq)) +
geom_boxplot()## Warning: Removed 46 rows containing non-finite values (stat_boxplot).
# look at diffs
beeAvgs$buzzDiff = beeAvgs$meanFreqRewarded - beeAvgs$meanFreqUnrewarded
hist(beeAvgs$buzzDiff)
Analysis with start/stop trials
GAMM4 shows time series
# start with gamm so I can show change by visit number
g01 = gamm4(freq ~ s(visitNum, by = trt2) + IT_mm + hive , random = ~(1|beeID), data = bees)
par(mfrow = c(2,2))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g01$gam) # Summary for paper ##
## Family: gaussian
## Link function: identity
##
## Formula:
## freq ~ s(visitNum, by = trt2) + IT_mm + hive
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 401.540 27.410 14.649 < 2e-16 ***
## IT_mm -16.875 6.467 -2.609 0.0091 **
## hive8 22.122 5.087 4.349 1.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum):trt2Rewarded 1st 50 Buzzes 4.550 4.550 28.068 <2e-16 ***
## s(visitNum):trt2Rewarded 2nd 50 Buzzes 1.332 1.332 1.267 0.359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0627
## lmer.REML = 59092 Scale est. = 1358.7 n = 5851summary(g01$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 59092
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8646 -0.4554 0.1523 0.6384 3.5375
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 551.509 23.484
## Xr.0 s(visitNum):trt2Rewarded 2nd 50 Buzzes 7.599 2.757
## Xr s(visitNum):trt2Rewarded 1st 50 Buzzes 657.613 25.644
## Residual 1358.685 36.860
## Number of obs: 5851, groups: beeID, 96; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 401.5401 27.4100 14.649
## XIT_mm -16.8750 6.4675 -2.609
## Xhive8 22.1223 5.0873 4.349
## Xs(visitNum):trt2Rewarded 1st 50 BuzzesFx1 -7.6413 5.6183 -1.360
## Xs(visitNum):trt2Rewarded 2nd 50 BuzzesFx1 0.8391 1.4259 0.588
##
## Correlation of Fixed Effects:
## X(Int) XIT_mm Xhive8 X(N15B
## XIT_mm -0.993
## Xhive8 0.097 -0.176
## X(N):2R150B 0.004 -0.004 0.006
## X(N):2R250B 0.002 0.001 0.037 0.001dev.off()## null device
## 1# save gamm plots
pdf(file.path(figDir, "Gamm_startStop_1st50.pdf"), width = 4, height = 3)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " frequency (Hz)"))
plot.gam(select = 1, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded first", side = 2, padj = -3.5)
abline(v=50, lty = 2)
dev.off()## null device
## 1pdf(file.path(figDir, "Gamm_startStop_2nd50.pdf"), width = 4, height = 3)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " frequency (Hz)"))
plot.gam(select = 2, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded second", side = 2, padj = -3.5)
dev.off()## null device
## 1LMER analysis
refref: todo – more model selection and diagnostics
– get p-vals, maybe?
mod1 <- lmer(freq ~ trt2 * treatment_rewarded + IT_mm + hive + (1|beeID), data = bees)
summary(mod1)## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 * treatment_rewarded + IT_mm + hive + (1 | beeID)
## Data: bees
##
## REML criterion at convergence: 59109.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8517 -0.4529 0.1408 0.6340 3.5575
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 540.1 23.24
## Residual 1367.5 36.98
## Number of obs: 5851, groups: beeID, 96
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 418.664 27.882
## trt2Rewarded 2nd 50 Buzzes -16.160 5.091
## treatment_rewarded T -12.449 1.301
## IT_mm -18.335 6.479
## hive8 23.248 5.074
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T 14.994 2.404
## t value
## (Intercept) 15.016
## trt2Rewarded 2nd 50 Buzzes -3.174
## treatment_rewarded T -9.566
## IT_mm -2.830
## hive8 4.582
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T 6.238
##
## Correlation of Fixed Effects:
## (Intr) t2R25B trtm_T IT_mm hive8
## trt2Rwr250B -0.227
## trtmnt_rwrT -0.046 0.176
## IT_mm -0.990 0.150 0.017
## hive8 0.122 -0.131 -0.024 -0.192
## tr2R250B:_T 0.012 -0.147 -0.542 0.002 0.041plot(mod1)
# diagnostics -- use REML = TRUE
m1 <- update(mod1, .~., REML =TRUE)
summary(m1) # summary for paper## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID) +
## trt2:treatment_rewarded
## Data: bees
##
## REML criterion at convergence: 59109.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8517 -0.4529 0.1408 0.6340 3.5575
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 540.1 23.24
## Residual 1367.5 36.98
## Number of obs: 5851, groups: beeID, 96
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 418.664 27.882
## trt2Rewarded 2nd 50 Buzzes -16.160 5.091
## treatment_rewarded T -12.449 1.301
## IT_mm -18.335 6.479
## hive8 23.248 5.074
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T 14.994 2.404
## t value
## (Intercept) 15.016
## trt2Rewarded 2nd 50 Buzzes -3.174
## treatment_rewarded T -9.566
## IT_mm -2.830
## hive8 4.582
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T 6.238
##
## Correlation of Fixed Effects:
## (Intr) t2R25B trtm_T IT_mm hive8
## trt2Rwr250B -0.227
## trtmnt_rwrT -0.046 0.176
## IT_mm -0.990 0.150 0.017
## hive8 0.122 -0.131 -0.024 -0.192
## tr2R250B:_T 0.012 -0.147 -0.542 0.002 0.041# pval for paper
m22 <- update(m1, .~. - trt2 : treatment_rewarded , REML = TRUE)
anova(m1, m22) # p-val for interaction## refitting model(s) with ML (instead of REML)## Data: bees
## Models:
## m22: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID)
## m1: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID) +
## m1: trt2:treatment_rewarded
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m22 7 59187 59234 -29587 59173
## m1 8 59150 59204 -29567 59134 38.864 1 4.544e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(m1)
qqnorm(ranef(m1)$beeID[[1]])
qqline(ranef(m1)$beeID[[1]])
sjp.lmer(m1, type = "re", sort = TRUE) # plot random effects to find any outliers## Plotting random effects...
sjp.lmer(m1, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.## Computing p-values via Wald-statistics approximation (treating t as Wald z).
Bootstrap CI’s for figure for paper
# set number of bootstrap replicates for models
nbootSims = 10
table(bees$hive) # more trials from hive 3##
## 7 8
## 3887 1964# using hive 7, since it's the one with the most data
# calculate an average IT for prediction
ITmean = mean(tapply(bees$IT_mm, INDEX = bees$beeID, FUN = function(x) x[1] ))
pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
IT_mm = ITmean,
hive = factor(7, levels = levels(bees$hive)),
beeID = 99999)
pframe$freq <- 0
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 10"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe## trt2 treatment_rewarded IT_mm hive beeID freq
## 1 Rewarded 1st 50 Buzzes F 4.267292 7 99999 0
## 2 Rewarded 2nd 50 Buzzes F 4.267292 7 99999 0
## 3 Rewarded 1st 50 Buzzes T 4.267292 7 99999 0
## 4 Rewarded 2nd 50 Buzzes T 4.267292 7 99999 0
## blo bhi predMean
## 1 332.3320 345.6363 340.4257
## 2 315.8471 328.5629 324.2655
## 3 322.0536 333.3766 327.9763
## 4 318.6930 331.4843 326.8099Make frequency plots for paper
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
facLevs = c("Rewarded 1st 50 Buzzes\n (buzz #51-100)",
"Rewarded 2nd 50 Buzzes\n (buzz #1-50)",
"Rewarded 1st 50 Buzzes\n (buzz #1-50)",
"Rewarded 2nd 50 Buzzes\n (buzz #50-100)"
)
pframe$t2 <- factor(facLevs, levels = facLevs[c(3,1,2,4)])
g00 <- ggplot(pframe, aes(x=t2, color = treatment_rewarded, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Treatment") +
theme(legend.position = c(.74,.8),
legend.background = element_rect(fill=alpha('gray90', 0.5)),
axis.text.x = element_text(angle = 25, hjust = 0.9)) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded?", begin =0.3, end = 0.8) +
geom_vline(aes(xintercept = 2.5))
g00
ggsave(plot = g00, filename = file.path(figDir, "StartStopPreds_freq2.pdf"), width =5, height = 3.5)
fac2Levs = c("51-100","1-50", "1-50","51-100")
pframe$buzzNums <- factor(fac2Levs, levels = fac2Levs[c(3,1)])
pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from = c(" F", " T"),
to = c("No", "Yes"))
g001 <- ggplot(pframe, aes(x=buzzNums, color = beesRewarded, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Sonication Number") +
theme(legend.position = c(.74,.8),
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded", begin =0.3, end = 0.8) +
facet_wrap(~trt2) +
ylim(c(315, 357))
g001
ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_freq3.pdf"), width =5, height = 3.5)
g0 <- ggplot(pframe, aes(x=trt2, color = treatment_rewarded, y=predMean))+
geom_point(position = position_dodge(width = 0.4))+
labs(y = "Sonication frequency (Hz)", x = "Treatment") +
theme(legend.position = c(.74,.8),
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1, position = position_dodge(width = 0.4)) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded", begin =0.3, end = 0.8) +
geom_vline(aes(xintercept = 1.5))
g0
ggsave(plot = g0, filename = file.path(figDir, "StartStopPreds_freq.pdf"), width =4, height = 3)Amplitude analysis for paper
mod0 <- lmer(log(amp_acc) ~ trt2 * treatment_rewarded + IT_mm + hive + (1|beeID), data = bees, REML = FALSE)
mod1 <- lmer(log(amp_acc) ~ trt2+ treatment_rewarded + IT_mm + hive + (1|beeID), data = bees, REML = FALSE)
BIC(mod1, mod0) # keep mod1## df BIC
## mod1 7 11347.39
## mod0 8 11354.21anova(mod0, mod1) #keep mod1## Data: bees
## Models:
## mod1: log(amp_acc) ~ trt2 + treatment_rewarded + IT_mm + hive + (1 |
## mod1: beeID)
## mod0: log(amp_acc) ~ trt2 * treatment_rewarded + IT_mm + hive + (1 |
## mod0: beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 7 11301 11347 -5643.3 11287
## mod0 8 11301 11354 -5642.4 11285 1.8507 1 0.1737mod2 <- update(mod1, .~. - trt2)
BIC(mod1, mod2) # keep mod2## df BIC
## mod1 7 11347.39
## mod2 6 11342.66anova(mod1, mod2) # disagrees with BIC## Data: bees
## Models:
## mod2: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
## mod1: log(amp_acc) ~ trt2 + treatment_rewarded + IT_mm + hive + (1 |
## mod1: beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod2 6 11303 11343 -5645.3 11291
## mod1 7 11301 11347 -5643.3 11287 3.9494 1 0.04689 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2 <- update(mod2, .~., REML = TRUE)
summary(m2) # summary for paper## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
## Data: bees
##
## REML criterion at convergence: 11310.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5886 -0.4828 0.1682 0.6921 2.4185
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.04746 0.2179
## Residual 0.39082 0.6252
## Number of obs: 5851, groups: beeID, 96
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.41964 0.26911 8.991
## treatment_rewarded T -0.12117 0.01804 -6.718
## IT_mm 0.35451 0.06346 5.586
## hive8 -0.25227 0.05043 -5.002
##
## Correlation of Fixed Effects:
## (Intr) trtm_T IT_mm
## trtmnt_rwrT -0.035
## IT_mm -0.992 0.004
## hive8 0.103 0.024 -0.182# pval for paper
m22 <- update(m2, .~. - treatment_rewarded , REML = TRUE)
anova(m2, m22) # p-val for interaction## refitting model(s) with ML (instead of REML)## Data: bees
## Models:
## m22: log(amp_acc) ~ IT_mm + hive + (1 | beeID)
## m2: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m22 5 11346 11379 -5667.8 11336
## m2 6 11303 11343 -5645.3 11291 44.929 1 2.043e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(m2)
qqnorm(ranef(m2)$beeID[[1]])
qqline(ranef(m2)$beeID[[1]])
sjp.lmer(m2, type = "re", sort = TRUE) # plot random effects to find any outliers## Plotting random effects...
sjp.lmer(m2, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects## Computing p-values via Wald-statistics approximation (treating t as Wald z).
Bootstrap Amp CI’s for figure for paper
# set number of bootstrap replicates for models
nbootSims = 100
table(bees$hive) # more trials from hive 3##
## 7 8
## 3887 1964# using hive 7, since it's the one with the most data
pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
IT_mm = ITmean,
hive = factor(7, levels = levels(bees$hive)),
beeID = 99999)
pframe$amp_acc <- 0
m00 <- update(mod0, .~., REML = TRUE)
pp <- exp(predict(m00, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m00, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 100"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe## trt2 treatment_rewarded IT_mm hive beeID amp_acc
## 1 Rewarded 1st 50 Buzzes F 4.267292 7 99999 0
## 2 Rewarded 2nd 50 Buzzes F 4.267292 7 99999 0
## 3 Rewarded 1st 50 Buzzes T 4.267292 7 99999 0
## 4 Rewarded 2nd 50 Buzzes T 4.267292 7 99999 0
## blo bhi predMean
## 1 49.00457 57.96004 53.73907
## 2 44.26685 52.01662 47.85855
## 3 42.39938 50.11380 46.54231
## 4 39.75753 48.64500 43.76588Make amplitude plots
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
facLevs = c("Rewarded 1st 50 Buzzes\n (buzz #51-100)",
"Rewarded 2nd 50 Buzzes\n (buzz #1-50)",
"Rewarded 1st 50 Buzzes\n (buzz #1-50)",
"Rewarded 2nd 50 Buzzes\n (buzz #50-100)"
)
pframe$t2 <- factor(facLevs, levels = facLevs[c(3,1,2,4)])
g00 <- ggplot(pframe, aes(x=treatment_rewarded, color = treatment_rewarded, y=predMean))+
geom_point()+
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Bee Rewarded") +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded?", begin =0.3, end = 0.8)
g00
ggsave(plot = g00, filename = file.path(figDir, "StartStopPreds_amp2.pdf"), width =5, height = 3.5)
fac2Levs = c("51-100","1-50", "1-50","51-100")
pframe$buzzNums <- factor(fac2Levs, levels = fac2Levs[c(3,1)])
pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from = c(" F", " T"),
to = c("No", "Yes"))
g001 <- ggplot(pframe, aes(x=buzzNums, color = beesRewarded, y=predMean))+
geom_point()+
labs(y =expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication Number") +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded", begin =0.3, end = 0.8) +
facet_wrap(~trt2)
g001
ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_amp3_mod0.pdf"), width =5, height = 3.5)
g001 <- ggplot(pframe, aes(x=beesRewarded, color = beesRewarded, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Sonication Number") +
theme(legend.position = c(.74,.8),
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded", begin =0.3, end = 0.8)
g001
bootstrap for a different model
# set number of bootstrap replicates for models
nbootSims = 5000
table(bees$hive) # more trials from hive 3##
## 7 8
## 3887 1964# using hive 7, since it's the one with the most data
pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
IT_mm = ITmean,
hive = factor(7, levels = levels(bees$hive)),
beeID = 99999)
pframe$amp_acc <- 0
m00 <- update(m2, .~., REML = TRUE)
pp <- exp(predict(m00, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m00, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 5000"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe## trt2 treatment_rewarded IT_mm hive beeID amp_acc
## 1 Rewarded 1st 50 Buzzes F 4.267292 7 99999 0
## 2 Rewarded 2nd 50 Buzzes F 4.267292 7 99999 0
## 3 Rewarded 1st 50 Buzzes T 4.267292 7 99999 0
## 4 Rewarded 2nd 50 Buzzes T 4.267292 7 99999 0
## blo bhi predMean
## 1 47.85513 54.45403 51.03044
## 2 47.85513 54.45403 51.03044
## 3 42.39319 48.27959 45.20699
## 4 42.39319 48.27959 45.20699pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from = c(" F", " T"),
to = c("No", "Yes"))
g001 <- ggplot(pframe, aes(x=beesRewarded, color = beesRewarded, y=predMean))+
geom_point()+
labs(y =expression ("Sonication amplitude "(m~s^{-2})), x = "Bee rewarded") +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray90', 0.5))) +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
scale_color_viridis(discrete = TRUE, name = "Bees Rewarded", begin =0.3, end = 0.8)
g001
ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_amp3_finalMod.pdf"), width =5, height = 3.5)Amplitude GAMM4
GAMM4 shows time series
# start with gamm so I can show change by visit number
g01 = gamm4(log(amp_acc) ~ s(visitNum, by = trt2) + IT_mm + hive , random = ~(1|beeID), data = bees)
par(mfrow = c(2,2))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, residuals = TRUE,pch = 20)
summary(g01$gam) # Summary for paper ##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum, by = trt2) + IT_mm + hive
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.40245 0.26419 9.094 < 2e-16 ***
## IT_mm 0.34419 0.06234 5.521 3.51e-08 ***
## hive8 -0.23957 0.04982 -4.808 1.56e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum):trt2Rewarded 1st 50 Buzzes 2.821 2.821 17.957 2.05e-10 ***
## s(visitNum):trt2Rewarded 2nd 50 Buzzes 2.150 2.150 3.144 0.0542 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0657
## lmer.REML = 11311 Scale est. = 0.39028 n = 5851summary(g01$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 11310.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5451 -0.4864 0.1642 0.6913 2.4609
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.04546 0.2132
## Xr.0 s(visitNum):trt2Rewarded 2nd 50 Buzzes 0.01600 0.1265
## Xr s(visitNum):trt2Rewarded 1st 50 Buzzes 0.02549 0.1596
## Residual 0.39028 0.6247
## Number of obs: 5851, groups: beeID, 96; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.402447 0.264191 9.094
## XIT_mm 0.344192 0.062339 5.521
## Xhive8 -0.239566 0.049822 -4.808
## Xs(visitNum):trt2Rewarded 1st 50 BuzzesFx1 0.026567 0.050437 0.527
## Xs(visitNum):trt2Rewarded 2nd 50 BuzzesFx1 -0.002741 0.046485 -0.059
##
## Correlation of Fixed Effects:
## X(Int) XIT_mm Xhive8 X(N15B
## XIT_mm -0.993
## Xhive8 0.106 -0.183
## X(N):2R150B 0.010 -0.010 0.012
## X(N):2R250B 0.000 0.002 0.028 0.001dev.off()## null device
## 1# save gamm plots
pdf(file.path(figDir, "Gamm_amp_startStop_1st50.pdf"), width = 4.5, height = 3.5)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " log amplitude ", (m~s^{-2})))
plot.gam(select = 1, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded first", side = 2, padj = -3.5)
abline(v=50, lty = 2)
dev.off()## null device
## 1pdf(file.path(figDir, "Gamm_amp_startStop_2nd50.pdf"), width = 4.5, height = 3.5)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " log amplitude ", (m~s^{-2})))
plot.gam(select = 2, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded second", side = 2, padj = -3.5)
dev.off()## null device
## 1